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1. INTRODUCTION 


Measuring the distance between geometric objects is an impor¬ 
tant problem in computer graphics, virtual reality, geomet ric mod- 
eling, computational geometry, CAD/CAM and robotics |Lin and| 
[Manocha 2003 1 . When objects are disjoint, the Euclidean distance 
between their closest points, also known as the separation dis¬ 
tance, is an obvious measure of distance. However, when objects 
overlap, the separation distance is undefined and a different mea¬ 
sure is needed to quantify the amount of interpenetration. Dif¬ 
ferent penetrati on measures have been intro duced, including pen¬ 
etration de pth ICameron and Culley 1986| , generalized pe netra¬ 
tion depth Zhang et al. 2007b|, pointwis epenetration depth | |Tang| 
|et al. 2009 , growth distance |Ong 1993] , and penetration volume 
IWeller and Zachmann 20091 . 


Penetration depth (PD) has been widely used by the computational 
geometry community. It is the distance that corresponds to the 
short est translation requir ed to separate two overlapping, rigid ob¬ 
jects IDobkin et al. 1993] . Many applications can benefit from PD 
computations. In physically-based animation, t he PD can be used 
to loc ate the point of application for impulses | |Guendelman et al.| 
|2003|, or used to find the time of contact in time-stepping meth- 
ods [Kim et al. 200^ . In constraint-based dynamics, when pene¬ 
tration is unavoidable due to numerical errors, PD can be used to 
roll back fro m an invalid state (penetration) to a valid one (non¬ 
penetration) |Redon 2004| . Moreover, PD can be used for post¬ 
stabilization such as the Baumgarte stabilization in rigid-body dy¬ 
namics to enforce non-penetration constrai nts. In virtual prot otyp- 
ing, the PD can be used to verify tolerances | |Kim et al. 2004a| . The 
length of a PD can be used to determine the appr opriate force feed - 
back in six-degree-of-freedom haptic systems [Kim et al. 2003 1 . 
Retraction-based motion planning algorithms can use PD to find 
a collision-free sample in t he narrow passage amongst obstacles 
[Zhang and Manocha 2008| , and PDs can be used to determ ine the 
existence of a path in planning scenarios [Zhang et al. 2008) . 

However, it is well-known that much more computation is gen¬ 
erally required to determine a PD than a separation distance. For 
general polyhedral objects, with a tota l of n faces, the computa¬ 
tion time is 0(n®) (Kim et al. 2004b| . This has motivated some 
researchers to find a more tractable approximation of the PD [Kim| 
|et al. 2004^ [Kim et al. 2002||Redon and Lin 2006} . Unfortunately, 
these methods are either too slow for interactive applications [Kim| 
|et al. 2004a[|Kim et al. 2002) or provide no guaranteed error bound 
[Redon and Lin 2066f. 

Main Contributions: We present a real-time algorithm to approx¬ 
imate the PD between arbitrary, polygon-soup models. Our algo¬ 
rithm actually determines an upper bound on the PD and we show 
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empirically that this is a tight bound on the exact value obtained 
from Minkowski sums. We also show how to obtain a set of local 
PD values from the PD solution. These local PDs characterize the 
amount of local overlap in each of several interpenetrating regions. 

Our algorithm is based on iterative and local optimization tech¬ 
niques. Given an in-collision configuration for an object, we find 
an initial collision-free configuration in configuration space. Then, 
we project this initial configuration on to a contact space using a 
variant of a continuous collision detection algorithm and construct 
a linear convex cone around the projected configuration. We then 
formulate the projection of the in-collision configuration onto this 
cone as a linear complementarity problem (LCP), which we solve 
using a Gauss-Seidel iteration. This runs until a locally optimal so¬ 
lution is obtained. 

Because ours is an iterative algorithm, finding a good initial config¬ 
uration is crucial. We therefore propose several search techniques 
to find the configuration from geometric properties such as the dis¬ 
tance between centroids, maximally clear configuration, sampling- 
based search or random search, as well as exploiting application- 
dependent information such as motion coherence. 

We have implemented our algorithm, PolyDepth, and benchmarked 
its performance in various complicated scenarios, such as random 
and pre-calculated configurations, and configurations governed by 
rigid-body dynamics. In all these scenarios, our algorithm achieves 
a highly interactive performance while providing a tight estimate 
of the PD value. 

Organization: The rest of the paper is organized as follows. In 
Sec.2, we briefly survey topics relevant to PD computations. We 
provide some preliminary information needed to understand our al¬ 
gorithm in Sec. 3, and give an overview of the algorithm. In Secs. 4 
and 6 we explain the two central techniques of our algorithm, out- 
and in-projection as well as local optimization techniques. In Sec. 
5 we explain how we find a good initial collision-free configura¬ 
tion. We present our experimental results, analyze the performance 
of PolyDepth in different benchmarks, and discuss implementation 
issues in Sec. 7. We conclude the paper in Sec. 8. 


2. PREVIOUS WORK 


|Kim et al. 2002| . Redon and Lin |Redon and Lin 200^ also pro¬ 
posed a CPU/GPU hybrid method of computing a lower bound on 
the P D, and this algorithm i s applicable to polygon-soup models. 
Lien (Lien 200^|Lien 2009) presented a sampling-based approach 
to PD computation base d on approximate M inkowski sum. More 
recently, Hachenberger (Hachenberger 2009) presented a method 
of obtaining an exact Minkowski sum based on convex decomposi¬ 
tion. However, these approaches are rather slow, and some require 
an explicit boundary representation of Minkowski sums that has to 
be re-evaluated whenever the orientation of an object is changed. 


A somewhat different route is to use distance fields that has the ad¬ 
vanta ge of being applicable to deformable models (Fisher and Lin| 
|2001| . Moreover, a GPU can be used to accelerate the computation 
of distance fields [Hoff et al. 2002[[Sud et al. 2006) . However, these 
approaches only provide a lower bound on a PD. There is also a 
variant of PD called a pointwise PD which is defined as the distance 
between the points of deepest interpenetration between t wo objects. 
Pointw ise PDs can be found using Hausdorff distance [Tang et ah] 
|2009| . However, a pointwise PD is merely a lower bound on the 
PD. A lower bound on a PD cannot be used to achieve separa¬ 
tion between overlapping objects, but an upper bound can. More 
recently, a penetra tion resolution tec hnique based on volume has 
been presented by [Allard et al. 2010[ . This method uses GPUs and 
can handle deformable objects. 


Generalized Penetration Depth: the constraints of some applica¬ 
tions mean that a translation is not sufficient to separate intersect¬ 
ing objects, and thus does not provide a useful measure of inter¬ 
penetration. In these cases, more complicated motions need to be 
considered. Zhang et al. proposed a generalized penetration depth, 
which is the minimal combination of translational and rotational 
motion need ed to separate two objects [Zhang et al. 2007b] |Zhang| 
|et al. 2007a| . More recently, kinematical geometr y has been used 
to compute generalized PDs [Nawratil et al. 200^. Non-Eucl i dean 
distan ce, such as growth distance [Ong and Gilbert 1996 |Ong| 
|1993| , can also be used as a measure of inter-penetration. Finally, 
som e researchers have investigat ed penetration volumes instead of 
PD [Weller and Zachmann 200^ , but the relation between this and 
the PD is questionable. 


There are several different types of PD algorithms for different 
model geometries and objectives. 

Convex Poly topes: a straightforward algorithm was presented to 
compute the PD b etween two convex polytopes by computing their 
Minkowski sums [Cameron and Culley 1986) . If the direction of 
motion is known, a minimal translational motion in this direc- 
tion can be found using a multi-resolution mesh hierarchy [Dobkin 
et al. 199^ . A randomized algorithm was presented by [Agarwal 
et al. 20001 . These algorithms provide exact solutions, but they 
are difficult to implement; in fact, no good implementations are 
known. However, various approximate algorithms have been de¬ 
velop ed based on upper and lower bounds on the PD [Cameron| 
|1997[ , expansion of p olyhedral approxim ations [Bergen 2001) , and 
dual-space expansion [Kim et al. 2004b| . The convex PD problem 
requires only O(n^) time in the worst case, where n is the total 
number of faces in the polytope. 

Non-convex Polyhedra: the computational complexity of PD for 
non-convex polyhedra is 0(n®), and thus no practical exact algo¬ 
rithms exist. A hierarchical refinement technique combined with 
GPU-assisted ray-shooting can provide an upper bound on the PD 


3. PRELIMINARIES 

We will now define the problem of penetration depth computation 
between general polygonal models and give an overview of our 
approach. 

3.1 Problem Formulation 

Suppose we have two objects 21 and 93 in R^. Without loss of gen¬ 
erality, we will assume that 21 is movable by translation and 93 is 
stationary, and both objects have the common global origin o. If 21 
and 93 are polyhedral objects and are interpe netrated, their penetr a- 
tion depth PD (21, 93) is formally defined as (Dobkin et al. 1993) : 

PD(2i, 93) = min{||d|| | interior(2i + d)n93 = 0, Vd G R^}. 

( 1 ) 

It is well known that the PD computation is closely related to the 
Minkowski sum. Formally, the Minkowski sums 21 0 93,21 0 —93 
between two compact sets, 21 and 93, are defined [Benson 1966| 
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ICameron 1997| as follows: 

21 © ^ = {a + b|a G 21, b G (2) 

21 0 -23 = {a - b|a G 21, b G 23}. (3) 

We can then reduce the problem of computing the PD between 
21 and 23 expressed by Eq. to the search for the minimum dis¬ 
tance betwee n o and the boundar y surface of their Minkowski sum, 
^(21©—23) [Dobkin etal. 1993) . However, since the interior of the 
polygon-soup models that we wish to consider may not be closed 
and thus be undefined, we need to define our PD problem in another 
way. 

Our definition of PD still follows the intuitive notion of using a min¬ 
imal translation to separate two overlapping models, even though 
the inside and outside of polygon-soup models is not properly de¬ 
fined. Thus, we define the PD to be the minimum distance from o 
to the boundary of the Minkowski sum 21 © —23: 

PD(2i, 23) = min{||d|| | interior(2l©—23)n{o+d} = 0,Vd G M^}. 

(4) 

Essentially, the boundary of the Minkowski sum constitutes the 
translational contact space between 21 and 23, given that 21 has a 
fixed orientation because it can only undergo translational motion. 

If 21 and 23 are polyhedra, then Eq[^is equivalent to Eq|^ 

3.2 Overview 

Although the exact PD can be computed from Minkowski sums, 
the explicit computation of Minkowski sums between complicated 
polygon-soup models is too slow for interactive applications. Look¬ 
ing at Eq. we see that finding the PD boils down to the search 
for the closest point to the origin on the contact space. Our algo¬ 
rithm approximates the contact space, which is the Minkowski sum 
boundary, only as far as is necessary to locate the closest point, 
refining its location until a locally optimal solution is obtained. 

The main computational components of our iterative algorithm are 
two projections: (a) a projection from an in-collision configuration 
on to the contact space (the in-projection); and (b) a projection from 
a collision-free or contact configuration on to the contact space (the 
out-projection). These two projections allow us to obtain a contact 
configuration from an in-collision configuration or from a collision- 
free configuration. 

The in-projection itself consists of two steps: (1) the construction 
of a local contact space (LCS) based on a linear complementar¬ 
ity problem (LCP) formulation; and (2) projection on to the LCS 
using a form of Gauss-Seidel iterative solver. The out-projection 


step is implemented using translational continuous collision detec¬ 
tion (CCD). Both out- and in-projections are explained in detail in 
Secs. 1^ an dm The overall flow of our iterative algorithm is as fol¬ 
lows (see Fig.[^: 

(1) Free-configuration selection: given two overlapping 
polygon-soup models 21 and 23, their common origin o 
should be inside the Minkowski sums 21 © —23. We select 
a collision-free configuration in the configuration space 
(Sec. 5). 

(2) Out-projection: then we perform out-projection from a source 
configuration = q^ to a target configuration q^ = o using 
CCD (Sec. 4). We call the configuration projected on to the 
contact space the current configuration q^. 

(3) In-projection: 

(a) We then find the contact features of q^, such as vertex/face 
(VE) or edge/edge (EE) contacts. From these features, we 
construct a local contact space (LCS) around q^, which is 
a linear convex cone in configuration space. 

(b) We perform in-projection from o to the LCS by formu¬ 
lating this problem as an LCP, which we solve using a 
form of Gauss-Seidel algorithm (Sec. 6.1). Thus this in- 
projected configuration becomes the new current projec¬ 
tion q^, and q^-i is set to the configuration obtained from 
the previous out-projection. 

(4) Sample classification: since q^ was obtained from the LCS, 
and not from the global contact space, q^ can be in-collision, 
in-contact or collision-free. We further classify the collision 
status of q^ by performing a static collision query on q^ as 
follows: 

(a) If q^ is an in-contact configuration, we compute the Eu¬ 
clidean distance from o to q^, and return it as the PD; i.e. 
PD(21, 23) = ||o — q^lj. The algorithm terminates. 

(b) If q^ is a collision-free configuration, q^ becomes the 
source configuration (q^ = Qi), and the origin becomes 
the target configuration (q^ = o). Then we go to step 2. 

(c) If q^ is an in-collision configuration, we obtain a proper 
contact configuration by setting q^ to a target configura¬ 
tion (q^ = q^), and the previous contact configuration 
q^-i becomes the source configuration (q® = Qi-i)- Then 
we go to step 2. 

(5) Iteration: steps 2-4 are repeated until the algorithm termi¬ 
nates. 
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(a) Local contact space and (b) Collision-free configura- (c) Out-projection from q-^ to (d) Construction of an LCS (e) In-projection on to the 

the origin tion q-^ o to obtain qo around qo LCS to obtain qi 

Fig. 2. Iterative optimization of a PD for a simple case of convex LCS. In this simple case, the PD algorithm terminates right after a single iteration of 
successive out- and in-projections. Then, PD = | |o — qi 11. 



LCS 


\ 



LCS 


\ 



qo/ 





LCS 


\ 


(a) Local contact space and the origin (b) Collision-free configuration q-^ (c) Out-projection from q-^ to o to (d) Construction of an LCS around 

obtain qo qo 


# • q^ • q^ • q^ 



(e) In-projection on to the LCS to ob- (f) Out-projection from qi to o to (g) LCS around q 2 (h) In-projection on to the LCS to ob¬ 
tain qi obtain q 2 tain qa 


Fig. 3. Iterative optimization of a PD for a slightly complicated case of convex LCS. In this case, PD requires two iterations since qi in (e) corresponds 
to separation. Then, PD = ||o — qa11. 


Note that this algorithm always terminates since it is a local op¬ 
timization on the LCS, and a locally optimal solution is obtained 
whenever q^ is in-conta ct. M ore discussion on termination condi¬ 
tions is provided in Sec. |6.3| In Figs.|^and[^ we illustrate our iter¬ 
ative algorithm for convex LCSs, and in Figs. and [TT] for more 
general cases. These examples illustrate the cases of translational 
configuration space for two dimensional polygonal objects. 


4. OUT-PROJECTION USING CONTINUOUS 
COLLISION DETECTION 

Our algorithm iteratively updates a sample configuration on the 
contact-space to find a locally optimal configuration. This up¬ 
date process requires a method of projecting the current in-contact 
or collision-free configuration on to another configuration on the 
contact-space. This out-projection is achieved by a variant of con¬ 
tinuous collision detection (CCD). 

4.1 Continuous Collision Detection 

Let 21 and 23 be two polygon-soup models in 3D, where 21 is mov¬ 
able by a translation M(t) and 23 is fixed. The source and target 
configurations of 21 are q^ and q^ at times t = 0 and t — 1, respec¬ 
tively. We also define 2l(t) = M(t)2l, and q(t) represents a con¬ 


figuration of 2l(t). Then the continuous collision detection (CCD) 
problem can be formulated as a search for the first time of contact 
(ToC) T, between [0,1], if it exists: 

T = min{t G [0,1] | 2l(t) fl 23 / 0}. (5) 

There are many met hods for CCD, but our choice is conservative 
advancement (CA) ||Lin 1993[ |Mirtich I996| |Zhang et al. 2006 1 
|Zhang et al. 2007c[|Tang et al. 2009| which is known to be fastest 

in practice. Like other CCD algorithms, CA takes the source q^ and 
target q^ configurations of an object 21, and computes the first time 
of contact when 21 linearly interpolates from q^ to q^ in the config¬ 
uration space. For a convex poly tope, CA computes a lower bound 
on T by repeatedly advancing 21 by At toward 23 until collision oc¬ 
curs. The value of At is chosen to correspond to a lower bound on 
the closest distance d(2l(t), 23) between 2l(t) and 23, and an up¬ 
per bound on the motion of 2i(t) projected on to the direction of 
(i(2l(t), 23) per second: 


At < 


d(2l(t),23) 


( 6 ) 


The first time of contact r is the sum of the time-steps At before 
collision; i.e. r — At. Whereas the CA algorithm applies to 
convex polytopes, we need to use the modified version of the C^A 
algorithm proposed by Tang et al. [Tang et al. 2009) for polygon- 
soup models. The C^A algorithm uses a bounding volume hierar- 


ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY. 


























PolyDepth: Real-Time Penetration Depth Computation using Iterative Contact-Space Projection 


5 


chy (B VH) based on swept sphere volumes (SSVs) |Larsen et aL| 
|2000) to control the depth of the recursive BVH traversal during it¬ 
erations of the algorithm, which reduces the computation time sig¬ 
nificantly. Since our PD problem is restricted to translational mo¬ 
tion, we can simplify C^A an d ma ke it faster. We present more 
details of this technique in Sec. |4.2| 

When we have foun d r, we can perform a static proximity query 
[Larsen et al. 2000| to find all the contact features between 21 (r) 
and 25, such as VF and EE contacts. These will be used to construct 
the boundary of the local contact space in Sec.[^ 

4.2 Translational Continuous Collision Detection 



Fig. 4. Minimal directional distance along v between two objects. 


as illustrated in Eig. Then the MDD between and is 
written A<b) = ^1|v||. We recursively compute the MDD 

between BV pairs or triangles pairs in the BVHs, and use Eq. □to 
obtain the ToC between the polygon-soup models 21 and 23. 

5. FINDING A COLLISION-FREE CONFIGURATION 

To obtain a non-trivial solution to Eq. the source configuration 
needs to be collision-free; otherwise, r is trivially zero. In the 
context of our PD computation, is unknown, whereas the target 
configuration q^, which is an in-collision configuration, is an input 
to the PD problem. We will introduce several methods of finding 
a collision-free source configuration q^. This is crucial to our PD 
computation, since our algorithm is a local optimization and starts 
from the contact configuration computed by CCD. 

5.1 Centroid Difference 

When no prior information is available about the interpenetration 
of the objects, the penetration direction can be estimated from the 
centroids of the objects. This direction can then be used to place 
the moving object in an interpenetration-free configuration, which 
can be expressed as follows: 


If 21 is moving with a constant translational velocity v, then we 
can call the shortest distance between 21 and 23 in the direction of 
V the minimal directional distance (MDD) (iv(2l, 23), as illustrated 
in Pig. The time at which 21 contacts 23 can be calculated as 
follows: 


dv(2l,23) 

l|v|| 


(7) 


If T < 1, then 21 and 23 will collide at r; otherwise, they are 
collision-free during the entire time-step. 

[Choi et al. 2006| presented a method of computing the MDD 
between convex polytopes based on Minkowski sums and ray¬ 
shooting. However, no algorithm for computing the MDD between 
polygon-soup models has been reported. We propose a simple 
method in which we construct the BVHs of polygon-soup models 
and recursively compute the MDD between pairs of nodes pairs in 
the BVHs; this is similar to the computation of Euclidean distance 
based on BVHs. Thus computing dv(2l, 23) boils down to com¬ 
puting the MDD between pairs of nodes in the BVHs (i.e. these 
nodes are bounding volumes (BVs) or triangles). We will explain 
how to compute the MDD between two triangles A^i and A<b ; and 
it should be apparent that a similar method can be used between 
BVs. 



Fig. 5. Translational CA for triangles. 


The ToC t' of two triangles under translational motion v can be 
found by repeated application of Eq. The motion bound is 
simply V • n, where n connects two closest points on the triangles. 


- o® 

q^ = q^ + (ra+ra) ||^a_o«|| ’ 

where q^ is the initial configuration of 21, and o® are the cen¬ 
troids of 21 and 23 respectively, and and r<B are the diameters of 
spheres enclosing each object. Eig.j^shows a collision-free config¬ 
uration determined from the centroid difference. 



Fig. 6. A collision-free configuration from the centroid difference. Ob¬ 
stacle 23 (yellow), initial in-collision configuration of 2t (magenta), and the 
collision-free configuration of 2t (red). 


5.2 Maximally Clear Configuration 

Although the centroid difference allows us to obtain a collision- 
free configuration, that configuration may be quite different from 
the optimal PD configuration for a complicated object with a lot 
of concavities or many holes. To obtain a collision-free configu¬ 
ration for objects of this sort, we find points in space that have 
maximal clearance. This preprocess allows objects to be positioned 
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(a) 



Fig. 7. Finding maximally clear configurations, (a) The given model, (b) AABB of the model, (c) Voxelize the AABB and compute the distance to the 
model from each grid position, (d) Remove the grid points corresponding to small distances, (e) Compare neighboring grid points and remove those with 
smaller distances by scanning along the x- and ^-directions, (f) Remove the grids on the boundary of AABB. (g) Compare neighboring grid points to remove 
those with smaller distances, and find maximally distant configurations in the x-, y- and ^-directions. 


without interpenetration. These points correspond to the points on 
the boundary of external Voronoi regions of the object, which are 
equidistant from at least two points on the surface of the model. 
These maximally clear configurations are appropriate candidates 
for collision-free configurations since they correspond to locally 
maximum likelihoods of a collision-free state (see Fig.[^. A similar 
strategy has bee n used in workspace sampling in motion planners 
ILatombe 199 1| . Since the construction of a generalized Voro noi 
diagram for a polygonal model is quite hard |Hoff et al. 1999) we 
propose a simple algorithm based on a voxelization of space (also 
see Fig.|^: 

(1) Compute the axis-aligned bounding box (AABB) of the static 
object, and voxelize the AABB with a grid. 

(2) Calculate distance to the object from each point in the grid. 

(3) Eliminate configurations with small distance values, since 
these nearly correspond to collisions. 

(4) Compare neighboring configurations and eliminate those with 
shorter distances; and find maximally distant configurations in 
one and two dimensions. 

(5) Remove any configurations remaining on the boundary of the 
AABB. 

(6) Again compare neighboring configurations and eliminate those 
with shorter distances, thus identifying the maximally distant 
configurations in full dimensions. 

Note that the above procedure only computes a subset of the points 
on the Voronoi boundary, which are those likely to have maximal 
clearance. 


(a) (b) (c) 

Fig. 8. Maximally clear configurations for the Cup, Grate, and 
Distorted-torus models. The small green spheres show the positions of 
the maximally clear configurations. 


5.3 Motion Coherence 

Applications of our PD algorithm are likely to exhibit motion co¬ 
herence: for instance, in rigid or articulated body dynamics, a se¬ 
ries of collision-free or contact configurations are calculated as a 
function of time; we can exploit the underlying motion coherence 
to find a good initial source configuration. A simple way of do¬ 
ing this is to cache a sequence of collision-free configurations and 
choose the one closest to a given in-collision configuration. More 
precisely, in our implementation, we store the last three in-contact 
configurations obtained from PD computations, and choose the one 
Qc that has the minimum translational difference from the input in¬ 
collision configuration o. Then, we change the orientation of Qc to 
be the same as that of o, and if this does not create any collision, 
we use Qc as a starting configuration for the iteration; otherwise, 
we slightly move Qc to the opposite direction toward o, and see 
if it causes collision again. If not, we use Qc as an initial configu¬ 
ration. Otherwise, we abandon the motion coherence strategy and 
switch to other methods (e.g. centroid difference). 

5.4 Random Configuration 

In a more general setting, in which we do not have any knowledge 
of q^, we can randomly sample in the free configuration space 
|Zhang et al. 2007 a| or simply locate q^ at infinity. 

5.5 Sampling-based Search 

The methods of finding a suitable collision-free configuration from 
which to begin out-projection that we have just described may still 
generate a configuration far from the optimal. A sampling-based 
search algorithm can be employed to refine an initial collision-free 
configuration. By performing collision detection on a series of con¬ 
figurations sampled on a line from the given collision configuration 
o to qo, as illustrated in Fig.M this sampling process terminates 
whenever a collision-free configuration is found. For collision de¬ 
tection, we use ILarsen et al. 2000| in our implementation. 

6. ITERATIVE OPTIMIZATION 

Any sampled configuration q in the contact space can be used as an 
estimate of the PD by computing its distance from the origin; i.e. 
PD(21, 23) = ||o — q||. When we have to deal with a highly non- 
convex object, finding a good candidate for q taxes the strategies 
suggested in Sec.[^ Thus we need to refine the sample to get a PD 
closer to the optimum, which we achieve by a walk in contact space. 
We use the result found by one of the techniques described in Sec. 
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Fig. 9. Sampling-based search for a collision-free configuration. The 

input collision configuration is o and Qq is an initial collision-free con¬ 
figuration. The first sample on the line from o to Qq is determined to be 

f 

in-collision, but the second sample is collision-free. 

[^as the start of this walk, and then build a local approximation to 
the contact space, and refine the configuration by in-projection. We 
repeat this process until a locally optimal solution is obtained. 

6.1 Contact-Space Localization and In-Projection 

The out-projection process explained in Sec. generates an in¬ 
contact configuration, which corresponds to a pair of contact fea¬ 
tures, one in each of 21 and 93. In general, this feature pair forms a 
vertex/face (VF), a face/vertex (FV) or an edge/edge (EE) primitive 
contact, from which we can construct a local contact space (LCS). 
Then, for the ith primitive contact, we determine the equation of a 
contact plane, j^q = where and respectively are the row 
vector of coefficients and the scalar bias of the zth plane equation. 

In some cases, this contact plane may degenerate to a lower¬ 
dimensional subset of configuration space, such as a line or a point. 
For example, contact between two collinear edges or a vertex/edge 
(VE) contact produces a contact line equation, and vertex/vertex 
(VV) contact produces a point. Our algorithm simply ignores these 
VV or collinear contacts, since projection on to a space of reduced 
dimensionality is unlikely to improve the PD. In the case of a VE 
contact, we generate several VE primitive contacts for all the faces 
that share the edge in the VE, and intersect them. The result is likely 
to be over-constrained s ince these contacts s hould be combined us¬ 
ing the union operator |Zhang et al. 2007^ . However, VE occurs 
rarely and the use of intersection simplifies implementation and im¬ 
proves performance. All other non-primitive contacts are also con¬ 
verted to several primitive contacts. 

Then, by stacking all the contact equations into a matrix, we can 
construct a system of linear equations: 

Jq = c. (8) 

We can then formulate in-projection as a minimization of the 
squared Euclidean distance between the origin o and a configu¬ 
ration q on the LCS, under the contact constraints: 

Minimize ||q||^ subject to Jq>c, (9) 

where Jq>c constrains q to lie either on the boundary of the LCS 
or outside it since the system of equations Jq = c represents the 
LCS. 

It is known [Redon et al. 2002| that Eq.j^is equivalent to a linear 
complementarity problem (LCP), formulated as a search for a value 


of A which satisfies the following complementarity condition: 

-|JJ^A + c > 0 

A > 0 (10) 

(-|JJ^A + c)^A = 0 

Then, q = | J^A. 

There are several me thods of solving an LCP such as Lemke’s and 
Dantzig’s algorithm [Cottle et al. 2009| ; however we choo se a pro¬ 
jected Gau ss-Seidel method for simplicity and stability [Jourdanj 
let al. 199^ . Then the LCP reduces to the search for a solution of 
the following linear system: 

JJ^A = 4c. (11) 

And the Gauss-Seidel iteration is: 

Afc+i = (D + L) ^(4c —UA/c), (12) 

where D + L + U = JJ^, and the matrices D, L and U respec¬ 

tively represent the diagonal, strictly lower triangular, and strictly 
upper triangular parts of the coefficient matrix JJ^, and k denotes 
the iteration number. When a component of A becomes negative 
during the iteration, we set it to zero: this is the projection step in 
the projected Gauss-Seidel (PGS) algorithm. 

A Gauss-Sei del iteration with an arbitrary symmetric positive def¬ 
inite matrix [Golub and Van Loan 199^ is known to converge. 
JJ^ is positive semidefinite since it is symmetric and JJ^x = 
(J^x)^J^x = ||J^x|p > 0 for X G We remove the 
redundant rows of J to make J full rank, and thereby promote 
JJ^ from semidefinite to positive definite since now x^JJ^x = 
(J^x)^J^x = ||J^x|p > 0 for all non-zero vectors x G 
Therefore, the Gauss-Seidel component of our algorithm always 
converges. 

We store the contact features as a list, sorted by distance from o 
to the corresponding contact plane. If a model is complicated, the 
number of contact features can become excessively high, and we 
select between only 10 and 30 contact feature pairs to reduce the 
complexity of solving Eq. This produces satisfactory results in 
practice. 

6.2 Iteration 

As described in the overview presented in Sec. |3.2| we iteratively 
optimize a sample to refine the penetration depth. Now we describe 
mor e details of this process, which relates to the five steps in Sec. 
[ 3:21 and Eig[g 

(a) An in-collision configuration (or the origin o) is given as input 
(Fig.[^a)). 

(b) Select a collision-free configura tion q^ using the technique de¬ 
scribed in Sec.[^(step 1 in Sec. |3.2| and Fig.p^b)). 

(c) Using the configurations q^ as source and o as target, perform 
out- proje ction to obtain a contact configuration qo (step 2 in 
Sec.|^and Fig.[^c)). 

(d) Construct an LC S around the contact configuration qo (step 
3(a) in Sec. |3.2| and Fig.[^d)). In this example, the LCS con¬ 
sists of only a single contact plane. 

(e) Perf orm in-projection on to the LCS to obtain qi (step 3(b) in 
Sec. |3.2| and Eig.p^e)). A proximity query is used to classify 
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(e) In-projection on to the LCS to ob- (f) Out-projection from qo to qi to (g) LCS around q 2 (h) In-projection on to the LCS to ob¬ 
tain qi obtain q 2 tain qa 


Fig. 10. Iterative optimization of a PD for a non-convex LCS. The PD is computed after two iterations since qi in (e) is in-collision. PD = | |o — qa 



(a) Collision at qa 


(b) Out-projection to obtain q 4 (c) In-projection to obtain qa 


Fig. 11. Iterative optimization of a PD for a more complicated LCS. Three iterations are required to obtain the PD since qa in (a) is in-collision. 
PD = ||o-q5||. 


the current configuration qi as contact, separation, or pene¬ 
tration. In this example, qi corresponds to penetration; so the 
next step is to find a valid in-contact configuration (step 4 in 
Sec.lT^. 

(f) Again perform out-projection to obtain q 2 (step 2 in Sec. |3.2| 
and Fig.p^f)). This requires a non-colliding source configura¬ 
tion as well as an in-collision target configuration. The source 
configuration (i.e. qo) was obtained from the previous out- 
projection, and the current configuration qi can be used as the 
target configuration. 


6.3 Algorithm Termination 

Our algorithm is based on an iterative optimization technique us¬ 
ing in- and out-projections. The out-projection ensures that our 
algorithm maintains an in-contact sample such that a proper up¬ 
per bound of PD is always obtained at any time of iteration. The 
in-projection ensures that our algorithm always reduces the PD 
value during iteration. These projections are iterated until a result 
from in-projection corresponds to an in-contact configuration so 
that no further reduction is possible. This way, our algorithm finds 
a locally-optimal solution for PD. 


(g) Construct the LCS around q 2 (step 3(a) in Sec. |3.2| and 
Fig.p^g)). This is a convex cone (the intersection of two con¬ 
tact planes), since q 2 consists of two contact feature pairs. 


(h) In-projection is performed again (step 3(b) in Sec. |3.2| and 
Fig.p^h)). The new configuration qa is in-cont act a nd the PD 
that will be returned is | jqa — o| | (step 4 in Sec.! 


3.2i. 


Since the complexity of the local contact space (LCS) is bounded 
above by 0(n®) where n is the number of polygons in polygon- 
soup models, the number of iterative projections in our algorithm 
is finite. Moreover, since our algorithm reduces a PD value at every 
iteration, no projection can be repeated for the same LCS. There¬ 
fore, our algorithm always terminates after a finite number of steps, 
which is typically 2 ~ 3 in practice. 


Fig. pT] shows a more complicated scenario. The first two steps ob¬ 
tain qi and q2, as before. But this time the proximity query at qa 
detects penetration (Fig.[TTJa)). Thus we perform out-projection to 
obtain q4 (Fig. [TTJb)), construct the LCS, and run in-projection 
again to get the final configuration qs (Fig.[TTJc)). 


6.4 Estimation of Local Penetration Depth 

When two objects intersect in multiple disjoint regions, applica¬ 
tions may require separate information about all the penetrations. 
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normals ni, n 2 

Fig. 12. Estimation of local PDs from the global PD. The red (53) and 
blue (2t) objects are placed in a penetrating configuration. The blue object is 
translated to the cyan one by the global PD d, and the two contact normals 
ni and n 2 are obtained at the contact. Projecting d on to ni and n 2 gives 
the two local PDs di and d 2 . 

i.e. local PDs instead of a single global PD. For instance, in rigid- 
body dynamics, a local PD, defined as the locally deepest penetrat¬ 
ing points o f two colliding objects |Guendelman et al. 2003 [[TSig] 
|et al. 2009), is required to compu te torques or to stabilize the simu- 
lation. [Guendelman et al. 2003 compute local PDs from distance 
fields, whereas [Tang et al. 2009 use two-sided Hausdorff distance. 

We propose an alternative method to derive local PDs from the 
global PD information that our algorithm computes. Our method of 
local PD estimation is based on the hypothesis that a pair of locally 
deepest penetrating points, one from each object, coincide when 
they are translated by the amount of the global PD. In other words, 
the locally deepest penetrating points are the ones that would be re¬ 
solved later than any other points; thus after the resolution by global 
PD, these points will just touch each other. Since PD allows only 
rigid translational motion, the amount of motion to resolve deep¬ 
est penetrating points is obtained by projecting the PD translation 
onto the penetrating direction of the penetrating points, which is 
identical to the normal of LCS corresponding to these points. 

More specifically, we compute the local PD of the zth interpen¬ 
etrating region as follows (also see Fig.p^: 

(1) Given two overlapping models 21 and 53, compute their global 
PD d = PD (21, 53) using the PD algorithm already presented. 

(2) Translate 21 by d to become 21(d), and compute the contact 
normal for each contact i between 21(d) and 53. 

(3) Project d on to each to get the local PD d^; i.e. d^ = (d • 

Note that our definition of local PD is not equivalent to p revious 
definition of [Guendelman et al. 2003 [ |Tang et al. 2009| , but it 
is computationally efficient. Our local PD algorithm is a simple 
by-product of our global PD algorithm, and does not require any 
costl y and memory-intensive precomputation of distance fields un¬ 
like [Guendelman et al. 2003), and no ex pensive Boolean intersec¬ 
tion is required unlike [Tang et al. 2009} . Also note that, according 
to our definition of local PD, even if d is not zero but if d is nor¬ 
mal to n, the local PD becomes zero. This is the case where the 
corresponding contact features are sliding. 


6.5 Computational Complexity of Algorithm 

Our algorithms consist of two main procedures: out- and in¬ 
projections. In our algorithm, the out-projection is implemented 
using translational CCD. The computational complexity of transla¬ 
tional CCD is similar to that of distance computation using bound¬ 
ing volume hierarchies (BVHs), since our translational CCD is 
equivalent to computing the m inimal directional d istance (MDD) 
between BVHs. According to [Larsen et al. 2000| , the cost func¬ 
tion T for distance calculation can be analyzed as T = x 
Cbv -\-NpXCp, where N^y is the number of bounding volume pair 
operations, C^y is the cost for pairwise distance computation for 
BVs, Np is the number of primitive pair operations, and Cp is the 
cost for pairwise distance computation for primitives. Translational 
CCD has a similar cost function, but C^y and Cp now correspond to 
the costs to compute the MDD for BVs and primitives, respectively. 
Moreover, we use an iterative method to calculate MDD, and each 
iteration requires Euclidean minimum distance calculation between 
a BV or triangle pair, which is a constant. In practice, the average 
number of iterations for MDD calculation is a small constant; for 
example, for two bunny models with 40K triangles each, as shown 
in Fig. it is 1.7. This means that the computational cost for 
translation CCD is only 1.7 times the cost for a minimum distance 
query (the function T above). 

The main computational cost for in-projection lies in solving 
Eq.[n] Since JJ^ is an n x n matrix where n is the number of 
contact features returned by out-projection, solving the linear sys¬ 
tem takes O(n^) time, but the Gauss-Seidel iterative solver would 
take O^Ngti^) time where Nq is the number of Gauss-Seidel iter¬ 
ations. Thus, the total computational complexity of our algorithm 
is (T + 0{NG'ri^))N where N is the total number of iterations for 
in- and out-projections. 



Fig. 13. Benchmarking models with triangle counts. Top row: Buddha 
(lOK), Bunnyl (IK), Bunny2 (40K), Cup (IK), Distorted-torus (1.33K), 
Dragon (174K), Golf-club (105K) Bottom row: Gratel (0.54K), Grate2 
(0.94K), Santa (152K), Spoon (1.34K), Torus (2K), and Torusknot (3K). 


7. RESULTS AND DISCUSSION 

We now present the results obtained by our PD computation algo¬ 
rithm in various scenarios, and discuss some of the implementation 
issues. 
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7.1 Implementation and Benchmarks 

We implemented our PD computation algorithnQ using Microsoft 
Visual C++ 8.0 under the Windows XP operating system. We tested 
the algorithm on a PC equipped with an AMD 2.22GHz CPU and 
2Gb of RAM. The benchmarking models that we used in these ex¬ 
periments are shown in Fig. The complexities of these models 
range from 0.54K to 174K triangles, and their topologies contain 
many holes and self-intersections (i.e. polygon-soups). We con¬ 
structed a number of benchmarking scenarios for these models, 
including random configurations, predefined sequences, and col¬ 
lisions within a rigid-body dynamics simulation. 

7.1.1 Random Configuration Scenario. We created 100 random 
configurations for several of the models, including the Torusknot, 
Buddha, Bunny2 and Dragon. We place a model at a fixed configu¬ 
ration, and translated a copy of the same model through a distance 
equal to half the size of its bounding volume, while changing its 
orientation, all at random creating many deeply penetrating con¬ 
figurations. Fig, shows an example in which two copies of the 
Torusknot model intersect. 






Fig. 14. Two Torusknot models in a random configuration. 


We use a centroid difference to find the initial collision-free config¬ 
uration. The performance of our algorithm is characterized in Table 
|T| These values are averaged for all the collision frames; the number 
of contacts means the number of contact feature pairs used to con¬ 
struct the LCS and the number of iterations is the total number of 
in- and out- projections for computing the PD. In Fig.[^ we also 
show the number of contact feature pairs, number of iterations, and 
computation time for two intersecting Torusknots in random con¬ 
figurations. 


Table I. Performance of our algorithm with random 
configurations. 


Model 

Time (msec) 

No. of Contacts 

No. of Iterations 

Torusknot 

5.53 

4.84 

2.81 

Buddha 

6.23 

5.04 

1.92 

Bunny2 

7.55 

8.25 

2.16 

Dragon 

12.76 

11.92 

2.29 


For the more topologically challenging models such as Grate 1, 
Grate2 and Distorted-torus, we use maximally clear configurations 
and sampling-based search scheme to find the initial collision-free 
configurations, because of their very concave geometries. Results 
are shown in Fig.p^ 

In order to assess the accuracy o f our PD algorithm, we used Lien’s 
open source library |Lien 200^ to compute near-exact Minkowski 


^The source codes are available for download at http://graphics. 
ewha.ac.kr/PolyDepth 


Fig. 15. Graphs of PD computation for two Torusknot models in 100 
random configurations. 

sums and PDs for the same models. Lien’s method is near-exact 
in the sense that it ignores low-dimensional boundaries of the 
Minkowski sums which have a volume close to zero. Fig.p^shows 
the PDs and the computational times for two Bunny 1 models and 
two Distorted-torus models using our method and Lien’s. 

The results are very close, even though our method is about 2000 
times faster than Lien’s. A quantitative comparison can be obtained 
by defining a metric for the relative error between approximate and 
exact PDs as follows: 


^PD ,relativel 


IIPD 


appoximate 


PD, 


0 — 55 ) 


, where w is the width of an object in a given direction n, defined 
as the distance between supporting planes normal to n. Since the 
PD is the minimum distance from the origin to the boundary of the 
Minkowski sums, maxu;(2i 0 —55) is an upper bound on the PD. 
However, since this metric depends on orientation, it can widely 
vary. An alternative metric which is unaffected by orientation is: 


D ,relative2 — 


approximate PD, 

2F(5l) 0 2f(55) 


, where v denotes the average magnitude of the vertices of 
an object, i.e. v = ^ ||vi||. If the objects are spheres, 

i 

^PD,reiativei = ^pD,reiative2, sincc the VS would bc thc radii 
of the spheres. If the PDs from Lien’s method are considered to be 
exact, then the errors in the results from our algorithm are listed in 
Tableini 
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Fig. 16. Using maximally clear configurations and sampling-based search to determine collision-free configurations for the Grate and Distorted- 
torus models. Obstacle 53 is yellow, and the input in-collision configuration of 21 is semitransparent magenta. The initial collision-free configuration of 21 is 
red, and the configuration of 21 translated by PD is shown in cyan. 





Fig. 17. Graphs of PD comparisons using our method and Lien’s 
method [Lien 200^ at random configurations. 

Table IT Average PD errors 


Model 

Mean error (%) 

Median error (%) 

Bunny1 

0.500 

0.165 

Distorted-Torus 

1.165 

0.066 



Fig. 18. Path for the Dragon model generated by a physics simulator. 



Fig. 19. Path of the Spoon model relative to Cup generated by a 
physics simulator. 

pairs, the number of iterations, and the computational time for the 
Dragon model. One Dragon is fixed in space and a copy falls under 
gravity, as shown in Fig. The centroid difference explained in 
Sec. |5.1| is used to find an initial collision-free configuration. 

Fig.[^shows scenes from the motion of the Spoon and Cup models 
along a predefined path. Collision-free configurations were found 
using the centroid difference (a^d) and using motion coherence 
(e^j). Fig. shows the results when maximally clear configura¬ 
tions were used for the same scene. Fig.[^a) shows the maximally 
clear configuration for the Cup. The centroid difference can gener¬ 
ate a poor collision-free configuration when the underling geometry 
and topology are complicated, as demonstrated in Fig.[^b). But 
in the same situations, the maximally clear configurations are very 
satisfactory, as demonstrated in Fig.[^e)^(j). The performance of 
our algorithm for predefined paths is summarized in Table [ni| 


Table III. PD performance on predefined paths 


Models 

Time (msec) 

No. of Contacts 

No. of Iterations 

Spoon and Cup 

0.96 

4.49 

1.02 

Buddha 

3.50 

5.07 

1.29 

Dragon 

5.84 

10.32 

1.45 


7.1.2 Predefined Path Scenario. Using Havok^^^ a rigid body 
dynamics simulator, we generated paths for the Dragon (Fig. p^ , 
and Spoon and Cup models (Fig.[^, and tested our PD algorithm 
along these paths. Fig. shows the number of contact feature 


^ http://www.havok.com 


7.1.3 Dynamics Scenario. We integrated our PD algorithm into 
rigid-body dynamics simulations involving the Torus, Bunny2, 
Golf-club and Santa models (see Fi g.[^. We used impulse-based 
dynamics | |Guendelman et al. 2003 1 to simulate rigi d-bod y dynam¬ 
ics based on the local PD method presented in Sec. |6.4| The local 
PDs are used to stablize the simulation and to resolve contacts and 
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Fig. 22. PD results for the Spoon and Cup models moving on predefined paths. Obstacles are shown in yellow, input in-collision configuration in semi¬ 
transparent magenta, initial collision-free configurations in red, contact configurations from in- and out-projections in green, and the configuration translated 
by the PD is shown in cyan. In (a) ~ (d), the centroid difference was used to find a collision-free configurations. In (e)~ (j), motion coherence was used. 



Fig. 23. Dynamic simulations results of pairs of Torus, Bunny2, Golf-club and Santa models using impulse-based dynamics IGuendelman et al.l 
[ 2 #^ . 


collisions. In this scenario, most collisions do not create deep pen¬ 
etrations since collisions are resolved immediately. As a result, the 
average PD computation is shorter, for example, than that required 
for a random configuration. We use maximally clear configurations, 
motion coherence and centroid differences to find initial collision- 
free configurations. A performance summary is given in Table [rv| 
Note that additional time is required for local PD computations. 


Table IV. Performance in the dynamics scenario. 


Models 

Global PD (msec) 

Local PD (msec) 

Total (msec) 

Torus 

3.13 

1.04 

4.17 

Bunny 2 

5.43 

1.78 

7.21 

Golf-club 

4.87 

1.67 

6.54 

Santa 

9.14 

2.90 

12.04 


7.2 Discussion 

There are some limitations to our algorithm. An object with high 
geometric and topological complexities may require a large num¬ 
ber of iterations to have a convergent solution. Furthermore, our 
method only approximates an upper bound, which can depend on 


the initial collision-free configuration, as Fig. [2^b) shows. Thus 
the quality of the initial collision-free configuration is critical. We 
have presented several methods of generating an initial collision- 
free configuration, and now summarize our experiences: 


(1) Maximally clear configurations are suitable for strongly non- 
convex objects of high genus, such as the Grate, Torus and Cup 
in Figs.f^andl^ since a collision-free configuration can be 


located inside a concavity. 


(2) Sampling-based search is useful for intricate or interlock¬ 
ing objects such as those shown in Fig. but it is time- 
consuming. 

(3) In a dynamics simulation, such as that shown in Fig.[^ motion 
coherence can be exploited to provide a good candidate for an 
initial configuration by using the configuration in the frame 
before a collision occurs. This can be achieved by caching. 

(4) The centroid difference is good for near-convex objects, unless 
they are deeply penetrating. In that case, we need to try several 
directions of back-off. 


(5) If the preprocessing cost of finding a maximally clear configu¬ 
rations is unaffordable, and no application-dependent informa¬ 
tion is available, random sampling can be employed. 
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(a) Number of Contact Features 



(b) Number of Iterations 



(c) Computation Time (msec) 


Nawratil et al.’s method | |Nawratil et al. 2009) iteratively calculates 
generalized penetration depth using kinematical geometry and con¬ 
strained optimization. However, their method has the following dif¬ 
ferences compared to ours: 

—Since the goal of their work is to compute a minimal rigid body 
displacement to separate overlapped objects by using a kinemat¬ 
ical mapping in their computational complexity is much 
higher than ours both in theory and in practice. It is unclear 
whether the application of this method to translational penetra¬ 
tion depth computation can perform at interactive rates like ours. 

—They also suggest two methods for finding an initial collision- 
free configuration, but this procedure takes a considerable 
amount of time to perform, e.g. a second for a model consisting 
of a few thousands triangles, which is too expensive for interac¬ 
tive applications. 

Allard et al .’s work is also related to penetration resolution |Allard| 
|et al. 20 T 0 I . However, their method is different from ours in the 
following sense: 

—Their work on penetration resolution is based on dynamics sim¬ 
ulation, whereas ours is based on purely geometric formulation, 
and they do not have any guarantees on collision resolution, 
whereas ours does. 

—Their work has no sense of optimality for collision resolution, 
whereas ours does. 

—Their work may suffer from the image resolution imposed by the 
use of GPUs, but ours does not. 

—Their work is suitable for deformable objects as well as rigid 
ones, whereas ours is suitable for rigid objects. 


Fig. 20. Graphs of PD computation using two colliding Dragon models 
moving on predefined paths. 



(a) (b) (c) 


Fig. 21. PD results using maximally clear configurations for Spoon 
and Cup models moving on predefined paths. 

In order to automate the selection of an initial configuration, we 
suggest following the above sequence and choose the resulting con¬ 
figuration closest to the input. However, none of these strategies 
guarantees a bound on the PD, since our algorithm uses local opti¬ 
mization; it may t ermina te far from the global optimum. However, 
the results of Sec. |7.1.l] show that these strategies work very well 
in practice, even for very challenging cases. 


8. CONCLUSIONS AND FUTURE WORK 


We have presented an interactive algorithm for computing the pen¬ 
etration depth of complicated polygon-soup models. Our method 
approximates the local contact space and iteratively performs in- 
and out-projections based on a Gauss-Seidel LCP solver and trans¬ 
lational continuous collision detection. We have proposed several 
schemes for selecting an initial collision-free configuration which 
utilizes motion coherence, centroid difference, and maximally clear 
configurations. We showed the effectiveness of our PD algorithms 
in various scenarios. We also presented a method of computing a 
local PD, and integrated it with a dynamics simulation. 


Future Work: We are interested in extending our algorithm to n- 
body PD problems, in which we need to separate multiple collid¬ 
ing bodies simultaneously. We would also like to extend our in- 
teractive algorithm to address the generalized PD proble m [Zhang | 
|et al. 2007b| jZhang et al. 2007a] Nawratil et al. 200^ in which 
both translational and rotational motions are possible. The key to 
making this problem amenable to our current framework is to find 
a way of linearizing the curved contact space. In addition, we are 
considering how our method might benefit other applications, such 
as haptic rendering and robot motion planning. 


7.3 Comparisons against Recent Approaches 

In this section, we make qualitative comparisons of our algorithm 
against re cent algorithms related to penetration de pth computation 
including | |Nawratil et al. 200^|Allard et al. 2010) . 
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